﻿// Copyright 2005, 2006 - Morten Nielsen (www.iter.dk)
//
// This file is part of SharpMap.
// SharpMap is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// SharpMap is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public License
// along with SharpMap; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 

// SOURCECODE IS MODIFIED FROM ANOTHER WORK AND IS ORIGINALLY BASED ON GeoTools.NET:
/*
 *  Copyright (C) 2002 Urban Science Applications, Inc. 
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

using System;
using System.Collections.Generic;
using System.Runtime.Serialization;
using Petra.Gis.Geometry.Support;
using Petra.Gis.ReferenceSystemServices;


namespace Petra.Gis.ReferenceSystemServices.Transformations
{
  /// <summary>
  /// Implements the Mercator projection.
  /// </summary>
  /// <remarks>
  /// <para>This map projection introduced in 1569 by Gerardus Mercator. It is often described as a cylindrical projection,
  /// but it must be derived mathematically. The meridians are equally spaced, parallel vertical lines, and the
  /// parallels of latitude are parallel, horizontal straight lines, spaced farther and farther apart as their distance
  /// from the Equator increases. This projection is widely used for navigation charts, because any straight line
  /// on a Mercator-projection map is a line of constant true bearing that enables a navigator to plot a straight-line
  /// course. It is less practical for world maps because the scale is distorted; areas farther away from the equator
  /// appear disproportionately large. On a Mercator projection, for example, the landmass of Greenland appears to be
  /// greater than that of the continent of South America; in actual area, Greenland is smaller than the Arabian Peninsula.
  /// </para>
  /// </remarks>
  [DataContract]
  public class MercatorProjection : ProjectionTransform
  {
    private double _falseEasting;
    private double _falseNorthing;
    private double lon_center;		//Center longitude (projection center)
    private double lat_origin;		//center latitude
    private double e, e2;			//eccentricity constants
    private double k0;				//small value m

    private double _semiMajor;
    //private double _semiMinor;

    private double _esinphi;

    /// <summary>
    /// xaml constr.
    /// </summary>
    public MercatorProjection() 
      :base()
    {
      //_className = "mercator";
    }




    /// <summary>
    /// Initializes the MercatorProjection object with the specified parameters to project points. 
    /// </summary>
    /// <param name="proj">projected coordinate system with Mercator projection</param>
    /// <remarks>
    /// </remarks>
    public MercatorProjection(ProjectedCoordinateSystem proj)
      : this(proj, false) { }




    /// <summary>
    /// Initializes the MercatorProjection object with the specified parameters.
    /// </summary>
    /// <param name="parameters">List of parameters to initialize the projection.</param>
    /// <param name="isInverse">Indicates whether the projection forward (meters to degrees or degrees to meters).</param>
    /// <remarks>
    /// <para>The parameters this projection expects are listed below.</para>
    /// <list type="table">
    /// <listheader><term>Items</term><description>Descriptions</description></listheader>
    /// <item><term>central_meridian</term><description>The longitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the longitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0).</description></item>
    /// <item><term>latitude_of_origin</term><description>The latitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the latitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0).</description></item>
    /// <item><term>scale_factor</term><description>The factor by which the map grid is reduced or enlarged during the projection process, defined by its value at the natural origin.</description></item>
    /// <item><term>false_easting</term><description>Since the natural origin may be at or near the centre of the projection and under normal coordinate circumstances would thus give rise to negative coordinates over parts of the mapped area, this origin is usually given false coordinates which are large enough to avoid this inconvenience. The False Easting, FE, is the easting value assigned to the abscissa (east).</description></item>
    /// <item><term>false_northing</term><description>Since the natural origin may be at or near the centre of the projection and under normal coordinate circumstances would thus give rise to negative coordinates over parts of the mapped area, this origin is usually given false coordinates which are large enough to avoid this inconvenience. The False Northing, FN, is the northing value assigned to the ordinate.</description></item>
    /// </list>
    /// </remarks>
    public MercatorProjection(ProjectedCoordinateSystem proj, bool isInverse)
      : base(proj, isInverse)
    {
      UpdateTransform();
    }






    ///// <summary>
    ///// Mercator projection axis list
    ///// </summary>
    //public override AxisList Axes
    //{
    //  get
    //  {
    //    // this is default Mercator axis list with correct orientation
    //    // this orientation must remain as defined because of projection equations
    //    return new AxisList()
    //    {
    //      new Axis("X", Axis.East),
    //      new Axis("Y", Axis.North)
    //    };
    //  }
    //}





    protected override bool parametersValid()
    {
      return Parameters.ContainsKey("central_meridian")
      && Parameters.ContainsKey("latitude_of_origin")
      && Parameters.ContainsKey("false_easting")
      && Parameters.ContainsKey("false_northing");

    }





    protected override void updateTransform()
    {
      //if (Projection.Projection.ClassName.ToLower() != "mercator_1sp" && Projection.Projection.ClassName.ToLower() != "mercator_2sp")
      //  throw new ArgumentException("invalid projected coordinate system, must have Mercator projection");

      //if (!(Projection.Unit is DistanceUnit))
      //  throw new ArgumentException(this.ToString() + " projection must have DistanceUnit type");

      //_projUnit = (DistanceUnit)Projection.Unit;

      Metadata = Projection.ProjectionDef.Metadata;

      _semiMajor = this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.SemiMajorAxis *
        this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.AxisUnit.MetersPerUnit;

      //_semiMinor = this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.SemiMinorAxis *
      //  this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.AxisUnit.MetersPerUnit;
      
      lon_center = Parameters["central_meridian"] * Projection.ProjectionDef.AngularParameterUnit.RadiansPerUnit;
      lat_origin = Parameters["latitude_of_origin"] * Projection.ProjectionDef.AngularParameterUnit.RadiansPerUnit;
      //ProjectionParameter scale_factor = _parameters["scale_factor"];
      _falseEasting = Parameters["false_easting"] * Projection.ProjectionDef.DistanceParameterUnit.MetersPerUnit;
      _falseNorthing = Parameters["false_northing"] * Projection.ProjectionDef.DistanceParameterUnit.MetersPerUnit;

      //double tmp = this._semiMinor / this._semiMajor;
      //e2 = 1 - tmp * tmp;
      //e = Math.Sqrt(e2);
      e = this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.Eccentricity;
      e2 = this.Projection.GeographicCoordinateSystem.Datum.Ellipsoid.e2;
      if (!Parameters.ContainsKey("scale_factor")) //This is a two standard parallel Mercator projection (2SP)
      {
        k0 = Math.Cos(lat_origin) / Math.Sqrt(1.0 - e2 * Math.Sin(lat_origin) * Math.Sin(lat_origin));
        //this.AuthorityCode = 9805;
        //this.Name = "Mercator_2SP";
      }
      else //This is a one standard parallel Mercator projection (1SP)
      {
        k0 = Parameters["scale_factor"];
        //this.Name = "Mercator_1SP";
      }
      //this.Authority = "EPSG";

    }






    /// <summary>
    /// Converts coordinates in decimal degrees to projected meters.
    /// </summary>
    /// <remarks>
    /// <para>The parameters this projection expects are listed below.</para>
    /// <list type="table">
    /// <listheader><term>Items</term><description>Descriptions</description></listheader>
    /// <item><term>longitude_of_natural_origin</term><description>The longitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the longitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0).  Sometimes known as ""central meridian""."</description></item>
    /// <item><term>latitude_of_natural_origin</term><description>The latitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the latitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0).</description></item>
    /// <item><term>scale_factor_at_natural_origin</term><description>The factor by which the map grid is reduced or enlarged during the projection process, defined by its value at the natural origin.</description></item>
    /// <item><term>false_easting</term><description>Since the natural origin may be at or near the centre of the projection and under normal coordinate circumstances would thus give rise to negative coordinates over parts of the mapped area, this origin is usually given false coordinates which are large enough to avoid this inconvenience. The False Easting, FE, is the easting value assigned to the abscissa (east).</description></item>
    /// <item><term>false_northing</term><description>Since the natural origin may be at or near the centre of the projection and under normal coordinate circumstances would thus give rise to negative coordinates over parts of the mapped area, this origin is usually given false coordinates which are large enough to avoid this inconvenience. The False Northing, FN, is the northing value assigned to the ordinate .</description></item>
    /// </list>
    /// </remarks>
    /// <param name="lonlat">The point in radians</param>
    /// <returns>Point in meters</returns>
    protected override void geographicToProjectedN(ICoordinate lonlat)
    {
      //if (double.IsNaN(lonlat.X) || double.IsNaN(lonlat.Y))
      //  return;
      //todo: invalid !!!! must use storage units def
      //double dLongitude = degrees2Radians(lonlat.X);
      //double dLatitude = degrees2Radians(lonlat.Y);

      //// convert to radians
      //lonlat.X *= _geogUnit.RadiansPerUnit;
      //lonlat.Y *= _geogUnit.RadiansPerUnit;

      /* Forward equations */
      if (Math.Abs(Math.Abs(lonlat.Y) - ProjMath.HALF_PI) <= ProjMath.EPSLN)
      {
        throw new ApplicationException("Transformation cannot be computed at the poles.");
      }
      else
      {
        _esinphi = e * Math.Sin(lonlat.Y);
        lonlat.X = _falseEasting + this._semiMajor * k0 * (lonlat.X - lon_center);
        lonlat.Y = _falseNorthing + this._semiMajor * k0 * Math.Log(Math.Tan(Math.PI * 0.25 + lonlat.Y * 0.5) * Math.Pow((1 - _esinphi) / (1 + _esinphi), e * 0.5));
      }

      //// convert to projcs units
      //lonlat.X /= _projUnit.MetersPerUnit;
      //lonlat.Y /= _projUnit.MetersPerUnit;

    }






    /// <summary>
    /// Converts coordinates in projected meters to decimal degrees.
    /// </summary>
    /// <param name="p">Point in meters</param>
    /// <returns>Transformed point in decimal degrees</returns>
    protected override void projectedToGeographicN(ICoordinate proj)
    {
      //double dLongitude = Double.NaN;
      //double dLatitude = Double.NaN;

      //// convert to meters
      //proj.X *= _projUnit.MetersPerUnit;
      //proj.Y *= _projUnit.MetersPerUnit;

      // todo: update parameters unit conversion
      //double dX = proj.X - this._falseEasting;
      //double dY = proj.Y - this._falseNorthing;
      proj.X -= this._falseEasting;
      proj.Y -= this._falseNorthing;


      /* Inverse equations
        -----------------*/
      double ts = Math.Exp(-proj.Y / (this._semiMajor * k0)); //t

      double chi = ProjMath.HALF_PI - 2 * Math.Atan(ts);
      double e4 = Math.Pow(e, 4);
      double e6 = Math.Pow(e, 6);
      double e8 = Math.Pow(e, 8);

      proj.lat = chi + (e2 * 0.5 + 5 * e4 / 24 + e6 / 12 + 13 * e8 / 360) * Math.Sin(2 * chi)
      + (7 * e4 / 48 + 29 * e6 / 240 + 811 * e8 / 11520) * Math.Sin(4 * chi) +
      +(7 * e6 / 120 + 81 * e8 / 1120) * Math.Sin(6 * chi) +
      +(4279 * e8 / 161280) * Math.Sin(8 * chi);

      //dLatitude = Phi2z(e,ts,out flag);
      /*if (flag != 0)
      {
        throw new ApplicationException();
      }*/
      proj.lon = (proj.X) / (this._semiMajor * k0) + lon_center;
      //dLongitude = AdjustLongitude(lon_center + dX/(this._semiMajor * k0));

      //proj.X = radians2Degrees(dLongitude);
      //proj.Y = radians2Degrees(dLatitude);

      //// convert to projcs units
      //proj.lon /= _geogUnit.RadiansPerUnit;
      //proj.lat /= _geogUnit.RadiansPerUnit;
    }




    /// <summary>
    /// Returns the inverse of this projection.
    /// </summary>
    /// <returns>IMathTransform that is the reverse of the current projection.</returns>
    public override SRSTransform Inverse()
    {
        return new MercatorProjection(Projection, !IsInverse);
    }




    public override Dimension TransformDimension
    {
      get { return Dimension.dim2D; }
    }




    public static string[] ProjectionNames
    {
      get { return new string[] { "mercator", "mercator_1sp", "mercator_2sp" }; }
    }




    //new public static MercatorProjection Create(ProjectedCoordinateSystem pcs)
    public static SRSTransform Create(ProjectedCoordinateSystem pcs, bool getAsInverse)
    {
      return new MercatorProjection(pcs, getAsInverse);
    }

  }
}
